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Theoretical uncertainties in the predictions of relativistic mean-field models are estimated using 
a chi-square minimization procedure that is implemented by studying the small oscillations around 
the chi-square minimum. By diagonalizing the matrix of second derivatives, one gains access to a 
wealth of information — in the form of powerful correlations — that would normally remain hidden. 
We illustrate the power of the covariance analysis by using two relativistic mean-field models: (a) 
the original linear Walecka model and (b) the accurately calibrated FSUGold parametrization. In 
addition to providing meaningful theoretical uncertainties for both model parameters and predicted 
observables, the covariance analysis establishes robust correlations between physical observables. 
In particular, we show that whereas the correlation coefficient between the slope of the symmetry 
energy and the neutron-skin thickness of Lead is indeed very large, a 1% measurement of the neutron 
radius of Lead may only be able to constrain the slope of the symmetry energy to about 30%. 
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I. INTRODUCTION 

The need to provide meaningful uncertainties in theoretical predictions of physical observables is a theme that is 
gaining significant momentum among the scientific community. Indeed, the search for a microscopic theory that both 
predicts and provides well-quantified theoretical uncertainties is one the founding pillars of the successful UNEDF 
collaboration [T]. Moreover, as recently articulated in an editorial published in the Physical Review A 2J, theoretical 
predictions submitted for publication are now expected to be accompanied by meaningful uncertainty estimates. The 
need for "theoretical error bars" becomes particularly critical whenever models calibrated in certain domain are used 
to extrapolate into uncharted regions. 

Although firmly rooted in QCD, computing both the nucleon-nucleon (NN) interaction and the properties of nuclei 
in terms of the underlying quark and gluon constituents remains a daunting task. Hence, rather than relying strictly 
on QCD, one uses the properties of QCD (such as chiral symmetry and relevant energy scales) as a guide to construct 
phcnomenological interactions using nucleons and mesons as the fundamental degrees of freedom. However, QCD has 
little to say about the strength of the underlying model parameters which must then be constrained from experimental 
data. For example, deuteron properties along with two-body scattering data are used to build a nucleon-nucleon 
interaction that may then be used (supplemented with a phenomcnological three-body force) to compute ab-initio 
the properties of light nuclei. Attempting ab-initio calculations of the properties of medium-to-heavy nuclei remains 
well beyond the scope of the most powerful computers to date. In this case one must bring to bear the full power 
of density functional theory (DFT) . Following the seminal work by Kohn and collaborators [3] , DFT shifts the focus 
from the complicated many-body wave-function to the much simple one-body density. Moreover, Kohn and Sham 
have shown how the one-body density may be obtained from a variational problem that reduces to the solution of a 
set of mean-field-like ( "Kohn-Sham") equations [3]. The form of the Kohn-Sham potential is in general reminiscent of 
the underlying (bare) NN potential. However, the constants that parametrize the Kohn-Sham potential are directly 
fitted to many-body observables (such as masses and charge radii) rather than two-body data. In this manner 
the complicated dynamics originating from exchange and correlation effects get implicitly encoded in the empirical 
constants. Yet regardless of whether the effective interaction is fitted to two-nucleon or to many-body data — the 
determination of the model parameters often relies on the optimization of a quality measure. 

In this contribution we focus on density functional theory and follow the standard protocol of determining the 
model parameters through a x 2 -minimization procedure. This procedure is implemented by: (a) selecting a set of 
accurately measured ground-state observables and (b) demanding that the differences between these observables and 
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the predictions of the model be minimized. Note that in the present framework a model consists of both a set of 
parameters and a x 2 -measure. In general, modifying the % 2 -measure (e.g., by adding observables) results in a change 
in the model parameters. Traditionally, once the % -minimum has been found one proceeds to validate the model 
against observables not included in the quality fit. Nuclear collective excitations are a potentially "safe" testing arena 
for the model as they represent the small oscillations around the variational ground state. But what happens when the 
model must be extrapolated to regions of large isospin imbalance and high density as in the interior of neutron stars? 
Clearly, without reliable theoretical uncertainties it is difficult to assess the predictions of the model. To remedy this 
situation we propose to study the small oscillations around the % -minimum — rather than the minimum itself. As we 
shall see, such a statistical analysis — inspired by the recent study reported in Ref. [5] — provides access to a wealth 
of information that remains hidden if one gets trapped in the x 2 -minimum. Among the critical questions that we 
will be able to answer is how fast does the x 2 -measure deteriorate as one moves away from the minimum. Should 
additional observables be added to the x 2 - m easure to better constrain the model? And if such observables are hard 
to determine are there others that may be more readily accessible and provide similar constraints? A particularly 
topical example that illustrates such a synergy is the correlation between the neutron-skin thickness and electric dipolc 
polarizability of neutron-rich nuclei [5 ,6. A detailed analysis of such correlation — which involves a systematic study 
of the isovector dipole response — is beyond the scope of this initial study and will become the subject of a forthcoming 
publication. Yet the study of correlations among observables sensitive to the poorly-determined density dependence 
of the symmetry energy will become a recurring theme throughout this contribution. 

The manuscript has been organized as follows. In Sec. [TT] we develop the necessary formalism to implement the 
correlation analysis. This section is divided in two parts: (a) a discussion on the structure of a class of relativistic 
mean-field models and (b) a relatively short — yet fairly complete — derivation of the statistical formalism required to 
perform the covariance analysis. In Sec. |III| two simple examples are used to illustrate the power of the formalism. 
This exercise culminates with the estimation of meaningful theoretical error bars and correlation coefficients. Our 



conclusions and outlook are presented in Sec. IV 



II. FORMALISM 



we 



In this section we develop the formalism required to implement the correlation analysis. First, in Sec. II A 
introduce a fairly general class of relativistic mean-field models that are rooted in effective-field-theory concepts, 
such as naturalness and power counting. Second, in Sec. |IIB] we present a self-contained derivation of the ideas and 
formulas required to implement the covariance analysis. 



Relativistic Mean-Field Models 



Relativistic mean-field models traditionally include an isodoublet nucleon field interacting via the exchange of 
two isoscalar mesons (a scalar </> and a vector one vector-isovector meson (6^), and the photon {A^) [THSj. The 
non-interacting Lagrangian density for such a model may be written as follows: 



1- ■ - • 1 2,2 



J% = <A (irfd»-M) $ + -d^d^ - -mict> 

- \V^V, V + \mlV»V, ■ V + \ml b" • b, \f^F, v , (1) 

where V^, b M „, and are the isoscalar, isovector, and electromagnetic field tensors, respectively. That is, 

V = dfjVv - d v V„ , (2a) 
b p „ = d^hv - d v b^ , (2b) 
Fp, = d„A v - d u Ap . (2c) 

The four constants M, m s , m v , and m p represent the nucleon and meson masses and may be treated (if wished) 
as empirical parameters. Often, however, m s is determined from an accurate calibration procedure. The interacting 
Lagrangian density has evolved significantly over the years and now incorporates a variety of meson self-interacting 
terms that are designed to improve the quality of the model. Following ideas developed in Ref. [3] we write the 
interacting Lagrangian density in the following form: 

J^int = $ [g a <t>- (sv^+l t ■ b M +|(i+r 3 )V) 7 M ] V - U(^V^b^) . (3) 
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In addition to the standard Yukawa interactions, the Lagrangian is supplemented with an effective potential 
U(<f>, Vfj,, b M ) consisting of non-linear meson interactions that serve to simulate the complicated dynamics that lies 
beyond the realm of the mean-field theory. Indeed, by fitting the various coupling constants directly to nuclear 
properties — rather than to two-nucleon data — the complicated dynamics originating from nucleon exchange, short- 
range effects, and many-body correlations gets implicitly encoded in a small number of parameters. For the purpose 
of the present discussion we introduce explicitly all non-linear terms up to fourth-order in the meson fields. That is, 

Ufa V\ b") = ^ + ^$ 4 - i (W^) 2 - A v (W V W) (B M • B^ 1 ) - 1 (B„ • B") * (4) 
+K <f>W tt W ,J '+K 1 <f>B fi ■ B^ + Xa^W^W + Xx^B^ ■ B» -k' v (w^W v ^ (V •B I/ )+... 

where the following definitions have been introduced: <& = g s (f), W ^ = g^V^, and B^ = g p h^. Given that the present 
analysis will be restricted to the study of uniform nuclear matter, terms proportional to the derivatives of the meson 
fields have not been included. As it stands, the relativistic model contains 14 undetermined parameters (1 meson 
mass, 3 Yukawa couplings, and 10 meson self-interaction terms). Note that if one incorporates the occasionally-used 
scalar-isovector ^-meson |10| lllj. then 9 additional parameters must be included to this order (1 Yukawa coupling, 
and 8 meson self-interaction terms). 

A model with 14 — or 23 — parameters goes significantly beyond the early relativistic models that were able to 
reproduce the saturation point of symmetric nuclear matter as well as various ground-state observables with only 
a handful of parameters (a single meson mass and three Yukawa couplings) J7J Q21 [13] . Although fairly successful, 
those early models suffered from a major drawback: an unrealistically large incompressibility coefficient. Such a 
problem was successfully solved by Boguta and Bodmer with the introduction of cubic and quartic scalar meson 
self-interactions [14]. Remarkably, using only these six parameters (m s , g s , g v , g p , k, X) it is possible to reproduce a 
host of ground-state properties of finite nuclei (both spherical and deformed) throughout the periodic table [T51 HE] ■ 
And by adding two additional parameters (£ and A v ) the success of the model can be extended to the realm of nuclear 
collective excitations and neutron-star properties [T7l420j . 

Given that the existent database of both laboratory and observational data appears to be accurately described 
by an 8-parameter model, is there any compelling reason to include 6 — or 15 — additional parameters? And if so, 
what criteria does one use to constrain these remaining parameters? A meaningful criterion used to construct an 
effective Lagrangian for nuclear-physics calculations has been proposed by Furnstahl, Scrot, and Tang based on the 
concept of "naive dimensional analysis" and "naturalness" [SUES]- The basic idea behind naturalness is that once 
the dimensionful meson fields (having units of mass) have been properly scaled using strong-interaction mass scales, 
the remaining dimensionless coefficients of the effective Lagrangian should all be "natural"; that is, neither too small 
nor too large [5TJ |33] . Such an approach is both useful and powerful as it allows an organizational scheme based on an 
expansion in powers of the meson fields. Terms in the effective Lagrangian with a large number of meson fields will 
then be suppressed by a large strong-interaction mass scale. In this regard the assumption of naturalness is essential 
as the suppression from the large mass scale should not be compensated by large, i.e., unnatural, coefficients. It 
was by invoking the concept of naturalness that we were able to truncate the effective potential U (<j>, V M , b M ) beyond 
quartic terms in the meson fields. 

Although we have justified the truncation of the effective Lagrangian invoking naturalness, we are not aware of an 
additional organizational principle that may be used a-priori to limit further the form of U(<p, V^,h^). This implies 
that all model parameters must be retained, as it is unnatural to set some coefficients arbitrarily to zero without 
a compelling symmetry argument |24j . In principle then, all model parameters must be retained and subsequently 
determined from a fit to empirical data. In practice, however, many successful theoretical models — such as NL3 [15U16] 
and FSUGold [TH] — arbitrarily set some of these parameters to zero. The "justification" behind these fairly ad-hoc 
procedure is that whereas the neglected terms are of the same order in a power-counting scheme, the full set of 
parameters is poorly determined by existing data, so ignoring a subset model parameters does not compromise the 
quality of the fit [9] 121). 

An important goal of the present work is to investigate correlations among the parameters of the model and 
whether additional physical observables could remove such correlations. To do so we follow the standard protocol of 
determining the model parameters through a x 2 minimization procedure. Traditionally, this procedure is implemented 
by selecting a set of accurately measured ground-state observables for a variety of nuclei and then demanding that 
the differences between the observables and the predictions of the model be minimized. Once this is done, the success 
of the model may be gauged by computing observables not included in the fit. However, it is often difficult to assess 
the uncertainty in the predictions of the model. To address this deficiency we propose to study the small oscillations 
around the minimum — rather than the minimum itself. Such a study — inspired by the recent statistical analysis 
presented in Ref. [5] — provides access to a wealth of information that, in turn, enables one to specify meaningful 
theoretical error bars as well as to explore correlations among model parameters and calculated observables. 
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Although the following discussion is framed in the context of an underlying x 2 -measure, our arguments are general 
as they merely rely on the existence of a (local) minimum (or an extremum). As in any small-oscillations problem, 
deviations of the x 2_m easure from its minimum value are controlled by a symmetric FxF matrix, where F represents 
the total number of model parameters. Being symmetric, such a matrix may be brought into a diagonal form by means 
of an orthogonal transformation. The outcome of such a diagonalization procedure is a set of F eigenvalues and F 
eigenvectors. When a point in parameter space is expanded in terms of these eigenvectors, the deviations of the x 2 ~ 
measure from its minimum value take the form of a system of F uncoupled harmonic oscillators — with the eigenvalues 
playing the role of the F spring constants. The spring constants may be "stiff" or "soft" depending on whether the 
curvature around the minimum is steep or shallow, respectively. As one explores the parameter landscape along a 
stiff direction — and thus along a particular linear combination of model parameters — a rapid worsening of the x 2 - 
measure ensues, suggesting that the fitting protocol is robust enough to constrain this particular linear combination. 
Conversely, no significant deterioration in the quality of the fit is observed as one moves along a soft direction. In 
this case the x 2 - mm imum is of little significance as scores of parameter sets (i.e., models) of nearly equal quality 
may be generated. This situation derives from the lack of certain critical observables in the x 2_m easure. As we shall 
see, the particular linear combination of model parameters defining the soft direction often provides enough hints 
to identify the missing observable(s). Moreover, through this sort of analysis one may establish correlations among 
observables that are particularly sensitive to such soft directions. This is important as certain observables may be 
easier to measure than others. A particular topical case is that of the neutron-skin thickness in 208 Pb and the electric 
dipole polarizability ® [35] . 



B. Linear Regression and Covariance Analysis 



As discussed earlier, relativistic models of nuclear structure are characterized by a number of model parameters, 
such as masses, Yukawa couplings, and non-linear meson coupling constants. Following the notation of Ref. [5J, we 
denote a point in such a parameter space by p= (pi, ■ ■ ■ ,Pf), where F is the total number of model parameters. In 
principle, each value of p represents a model. In practice, of course, one is ordinarily interested in the "best model" 
as defined by a quality measure. To do so, the model parameters are often calibrated to a well-determined set of 
ground-state properties of finite nuclei (such as masses and charge radii) that is supplemented by a few bulk properties 
of infinite nuclear matter (such as the binding energy, incompressibility coefficient, and symmetry energy at saturation 
density). Once the model parameters and the group of observables have been selected, the optimal parameter set is 
determined via a least-squares fit to the following x 2 quality measure: 

X (P) = >. 44 ■ (5) 




Here N (often much larger than F) denotes the total number of selected observables whereas "th" and "exp" stand 
for the theoretical prediction and experimental measurement, respectively. Further, every observable is weighted by 
a factor of (AO n )~ 1 that is (customarily) associated with the accuracy of the measurement. 

We assume that — through a numerical procedure that is not of particular relevance to this work — an accurately- 
calibrated model p has been found. This implies that all first derivatives of \ 2 vanish at p . That is, 



dx 2 {p) 



dpi 



=P0 



= <9 lX 2 ( Po ) = (fori=l,...,F). (6) 



The existence of the minimum (as opposed to a maximum or saddle point) also implies that a particular set of F 
second derivatives (to be defined shortly) must all be positive. Approaches based on a least-squares fit to a x 2 -measure 
often culminate with the identification of the optimal parametrization po. The predictive power of the model may 
then be appraised by computing observables that were not included in the fitting protocol. Less often, however, 
least-squares-fit approaches are used to evaluate the "uniqueness" of the model. In other words, how fast does the 
X 2 -measure deteriorate as one moves away from p ? Clearly, if the minimum is relatively flat (at least along one 
direction), then there will be little (or no) deterioration in the quality of the fit. Through a statistical analysis, we 
will be able to obtain a physically reasonable domain of parameters. We implement such analysis by studying the 
small oscillations around the x -minimum. As a bonus, we will be able to uncover correlations among observables 
and attach meaningful theoretical error bars to the theoretical predictions [5J. To start, we expand the ^-measure 
around the optimal po model. That is, 

! F 

X 2 (P) = X 2 (Po) + 2 (P ~ P°Mp " Po)A9 3 X 2 (Po) + ■■■ (7) 
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For convenience, we quantify the departure from the minimum by defining scaled, dimensionless variables 

= (P ~ Po)» 
(Po)i 



(8) 



In terms of these scaled variable, the quadratic deviations of the x 2 - m easure from its minimum value take the following 
compact form: 

X 2 (P) - X 2 (Po) = Ax 2 (x) = x T M x , (9) 

where x is a column vector of dimension F, x T is the corresponding transpose (row) vector, and M. is the symmetric 
Fx F matrix of second derivatives defined by 



Ma = 2 [g^: ) _ o = ^(poMpo W^(po) ■ (10) 

Being symmetric, the matrix M. can be brought to a diagonal form by means of an orthogonal (change-of-basis) 
transformation. Denoting by A the orthogonal matrix whose columns are composed of the normalized eigenvectors 
and by 2? = diag(Ai, . . . , \p ) the diagonal matrix of eigenvalues, the following relation holds true: M — AVA T . By 
inserting this relation into Eq. ([9]), we obtain the following simple and illuminating expression: 

a x 2 (x) = ^{Af>A T y = en = E ■ (ii) 

i=l 

Here the vector £ = A T x represents a point in parameter space expressed, not in terms of the original model parameters 
(ffsi 9vi - • •) but rather, in terms of the new ( "rotated") basis. As previously advertised, the deviations of the x 2 -measure 
from its minimum value have been parametrized in terms of F uncoupled harmonic oscillators — with the eigenvalues 
playing the role of the spring constants. In this way, each eigenvalue controls the deterioration in the quality of the 
fit as one moves along a direction defined by its corresponding eigenvector. A "soft" direction — characterized by a 
small eigenvalue and thus little deterioration in the \ 2 measure — involves a particular linear combination of model 
parameters that is poorly constrained by the choice of observables included in the least-squares fit. By isolating such 
linear combination(s) one can identify what kind of observables (e.g., isovector observables) should be added to the 
X 2 -measure to better constrain the theoretical model. Moreover, one may also explore correlations among various 
observables (e.g., neutron-skin thickness and dipole polarizability) thereby facilitating the experimental extraction of 
some of these critical observables. This could be done by either refining existing experimental measurements or by 
designing brand new ones. 

A concept of fundamental importance to the correlation analysis is the covariance between two observables A and 
B, denoted by cov(A, B) [35]. Assuming that (xW, . . . , x^ M )) represent M points (or models) in the neighborhood of 
the optimal model x' ^ =0, the covariance between A and B is defined as: 

M 

cov(A,B) = - E [(> n) ^ (A)) (S (m) - <£))] = (AB) (A)(B) , (12) 

m— 1 

where A^ = /l(x( m )) and "()" denotes a statistical average. From the above definition the correlation coefficient — 
often called the Pearson product-moment correlation coefficient — now follows: 

P (A )B )= ; ov{A ' B) , as) 

where the variance of A is simply given by var(A) = cov(A, A) . Note that two observables are said to be fully correlated 
if p(A, B) = 1, fully anti-correlated if p(A, B) = —l, and uncorrelated if p(A, B) = 0. If one expands the deviation of 
both observables from their average value, then cov(A, B) may be written as 



M 

(m) (m) 



M ^ Xl ' X i 

m— 1 



dB^ = y^dA c dB_ 

dxj ^ dxi 13 dx,j 
l j 



where both derivatives are evaluated at the minimum (x(°)=0) and the covariance matrix C'ij has been introduced 26 . 
In order to continue, it is critical to decide how should the M points be generated. A particularly convenient choice 



G 



is to assume that these M points (or models) are distributed according to the quality measure x 2 ( x )- That is, we 
assume a probability distribution </>(x) given by 



(x) = exp 



;A X 2 (x) 



= exp ( — -x M. x 



The covariance matrix may then be written as follows: 

_ J XiXj4>(x)dx _ 1 
lJ ~ J 0(x)dx ~ Z(0) 

where we have defined the "partition" function Z(3) as 



d 2 Z{3) 



dJidJj 



j=o 



Z(J) 



(x)e J ' x dx = / exp ( ~^x T M x + J ■ x ) dx 



The above gaussian integrals may be readily evaluated by completing the square. One obtains 

Z(J) = Z(0)cxp (-J* A- 1 j\ = Z(Q)e w{3) . 



(15) 



(16) 



(17) 



(18) 



Hence, under the assumption that the model parameters are generated according to the x 2 -measure, the covariance 
matrix becomes equal to the inverse of the matrix of second derivatives of x 2 . That is, 



C - — 
13 ~ Z(0) 



d 2 Z(3) 
dJidJj 



j=o 



dJ,d,L 1 h 



(19) 



Finally then, we arrive at a form for the covariance of two observables that is both simple and easy to compute: 



cov(A,B) = ]T 

i,3 



dxi 13 dxj 



1 dB_ 



(20) 



The last term in the previous expression is particularly illuminating. Consider, for example, the case of a very soft 
direction in the x 2 - m easure, namely, an eigenvector of M (say £j) with a very small eigenvalue (say A,^ 1 3> 1). Such a 
situation routinely emerges in RMF models whenever two or more isovector parameters are included in the Lagrangian 
density but only masses and charge — not neutron — radii are used to define the x 2_m easure. Having identified a soft 
direction, one could then search for an observable A (e.g., the neutron-skin thickness in 208 Pb) that is particularly 
sensitive to such a soft direction (as indicated by dA/d^i 1). Adding such an observable to the x 2 -measure will stiffen 
the formerly soft direction, thereby improving the predictive power of the model. Moreover, if A is difficult to measure, 
one could search for alternative observables that are strongly correlated to A. Although some of these notions have 
been heuristically implemented for some time, the statistical analysis discussed here provides a quantitative measure 
of the correlation between observables [S] . 



III. RESULTS 



In this section we provide two simple examples that illustrate the ideas presented in the previous sections. Here 
terms such as "unique" and "predictive" will be used to characterize a model. We regard a model as being unique 
if all the eigenvalues of Ai are large (i.e., A, ^ 1 for all i). A model is predictive if it can successfully account for 
physical observables not included in the x 2_m easure. Note that a model has been defined here as consisting of both 
an underlying Lagrangian density (or effective interaction) and a set of physical observables defining the x 2_m easure. 



A. Example 1: Linear Walecka Model 



We start this section by discussing the linear Walecka model as an example of a model that is unique but not 
predictive. The Lagrangian density for this case is simple as it only contains two coupling constants [71 112). That is, 



(21) 
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The Walecka model is perhaps the simplest model that can account at the mean-held level- for the saturation 
of symmetric nuclear matter. Indeed, it is the saturation density and the energy per nucleon at saturation that 
are typically used to calibrate the two parameters of the model. To make this simple model slightly less trivial we 
determine the two parameters of the model by minimizing a quality measure x 2 defined in terms of three "observables" : 
(i) the saturation density po, (ii) the energy per nucleon at saturation Sq, and (iii) the effective Dirac mass Mq. Central 
values and uncertainties for these three quantities are given as follows: 

po = (0.155 ± 0.01) fm~ 3 , (22a) 
so = (-16 ± 1) MeV , (22b) 
Mq = (0.6 ± 0.1) M . (22c) 

Using standard numerical techniques, a minimum value for the x 2 -measure f Xo — 0.34145 is obtained at 

gl = 93.62647 , (23a) 
gl = 180.48347 . (23b) 

Having computed the minimum value of the quality measure, we now examine its behavior arou nd the minimum. 



This is implemented by diagonalizing the symmetric matrix of second derivatives M. [see Eq. (10)]. The outcome of 
such a diagonalization procedure is the diagonal matrix of eigenvalues T> and the orthogonal matrix of normalized 
eigenvectors A. That is, 

V = diag(Ai,A 2 ) = diag(7.4399xl0 4 ,8.3195xl0 1 ) , (24a) 

2 _( cos 6 sine A _ / 0.74691 0.66492 A , . 

A ~ \ -sin9 cos 9 ) ~ \ -0.66492 0.74691 ) ' ^ > 

It is evident that both eigenvalues are very large. This indicates that both directions in parameter space are stiff 
and consequently the quality measure (Ax 2 = Ai£ 2 + A 2 £ 2 ) will deteriorate rapidly as one moves away from the x 2 - 
minimum. Note that Ai is significantly larger than A 2 ; this is to be expected. When probing the parameter landscape 
along the first direction (i.e., £ 2 =0) the scalar and vector coupling constants move "out-of -phase" (see the first column 
of the matrix A). For example, the scalar attraction would get larger at the same time that the vector repulsion 
would get smaller. This would yield a significant increase in the binding energy per particle and consequently a drastic 
deterioration in the % 2 -measure. Recall that large and cancelling scalar and vector potentials are the hallmark of 
relativistic mean-field models. 

To quantify the extent by which the linear Walecka model is unique, we now proceed to compute the variance in 



the coupling constants using Eq. ( 20 ) . We obtain 



(j 2 = = A^ 1 cos 2 9 + A^ 1 sin 2 6 = 5.3217 x 10~ 3 , (25a) 

cr 2 = (A^ 1 ) 22 = sin 2 9 + A^ 1 cos 2 9 = 6.7116 x 10~ 3 . (25b) 

In turn, this translates into the following uncertainties in the optimal values of the coupling constants: 

g 2 s = 93.62647 (1 ± a s ) = 93.62647 ± 6.83008 , (26a) 
gl = 180.48347 (1 ± cr v ) = 180.48347 ± 14.78596 . (26b) 

We conclude that the uncertainties in the model parameters — and thus in most of the predictions of the model — are 
of the order of 5-to-10 percent. In principle, the model uncertainties could be reduced by refining the experimental 



database [see Eq. (22 1]. The great merit of the present statistical approach is that one may systematically explore the 



extent by which the experimental measurement must be refined in order to achieve the desired theoretical accuracy. 



Note that the theoretical uncertainties are dominated by the smallest eigenvalue of M. [see Eq. (25)]. Thus, assessing 



the uniqueness of the model by varying each model-parameter individually (e.g, first gl and then gl) is misleading 
and ill advised. It is misleading because in doing so the quality measure will in general be dominated by the largest 



eigenvalue [see Eq. (11)]. Yet it is the lowest eigenvalue that determines the uniqueness of the model. 



Carrying out the covariance analysis further, we now proceed to compute correlation coefficients between model 



parameters and observables [see Eqs. (13) and (20)]. In estimating uncertainties in the model parameters one con- 



centrates on the diagonal elements of the (inverse) matrix of second derivatives [see Eq. (25)]. Information on the 



correlation between model parameters is, however, stored in the off-diagonal elements. For example, the correlation 
coefficient between g^ and g^ is given by 



P(fs.5v) = 



= 0.9977 . 



M- 1 ) [M- 1 ' 

11 V / 22 



(27) 



The strong (positive) correlation between gl and is easily understood. Given that configurations in parameter space 
are distributed according to the x 2 -measure, model parameters in which g^ and gl move out-of-phase are strongly 
suppressed, as they are controlled by the largest eigenvalue \\. As a result, an overwhelming number of configurations 
are generated with g^ and g^ moving in phase, thereby leading to a large positive correlation. Correlation coefficients 
between various isoscalar observables have been tabulated in Table Q] Given that the correlation coefficients are 



sensitive to the first derivatives of the observables along all (eigen) directions [see Eq. (20 1], we have listed them for 
completeness in Table [TTJ We observe that all observables display a much larger sensitivity to the stiff direction than 
to the soft one. This could (and does) lead to sensitive cancellations since the large derivatives compensate for the 
small value of A^f 1 . Indeed, the correlation between the saturation density and the binding energy at saturation is 
very small. On the other hand, the incompressibility coefficient appears to be strongly correlated to the binding 
energy. This behavior is also displayed in graphical form in Fig. [T] where predictions for the various observables were 
generated with model parameters distributed according to </>(x) [see Eq. (15)]. Note that the covariance ellipsoids in 
Fig. [TJwere generated by selecting those model parameters that satisfy A;p~< 1. 





eo 


Po 


K 


M * 


£(> 


+1.0000 


-0.0036 


-0.9998 


+0.8867 


Po 


-0.0036 


+1.0000 


-0.0174 


+0.4591 


K 


-0.9998 


-0.0174 


+1.0000 


+0.8962 


M * 


+0.8867 


+0.4591 


+0.8962 


+1.0000 



TABLE I: Correlation coefficients between isoscalar observables in the linear Walecka model. 




520 540 560 580 0.14 0.15 0.16 0.17 



K (MeV) Po (fm - 3 ) 

FIG. 1: (Color online) Predictions from the linear Walecka models for the saturation density, binding energy, and incompress- 
ibility coefficient at saturation. Modei parameters were generated according to the distribution exp(— Ax 2 /2). Both of the 
covariance ellipsoids were generated by limiting the models to the region A% 2 < 1. 



Based on the previous statistical analysis it appears that the linear Walecka model is unique (at least at the 5 — 10% 
level). But is the linear Walecka model predictive? To test the predictability of the model we focus on two physical 
observables that were not included in the x 2 -uieasure, namely, the incompressibility coefficient Kq and the symmetry 
energy J. We obtain — with properly computed theoretical errors — the following results: 



K = (552.537+ 29.655) MeV , 
J = (19.775 ± 0.683) MeV . 



(28a) 
(28b) 
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de 


-1.6698X10 1 


-1.1716xl0 _i 


dpo 


3.6619 


-5.7333 xlO -1 


dK 


1.4261 xlO 1 


1.1055 xlO" 1 


dMS 


-3.2349 


-8.8817 XlO - * 



TABLE II: First derivatives of the scaled observ able s (i.e., observable scaled to its value at the x 2 -minimum) as a function of 



£i and £2 evaluated at the \ -minimum; see Eq. (201 



Both predictions, even after theoretical errors have been incorporated, differ significantly from the presently acceptable 
values of K ss (240 ± 20) MeV and J w (32 ± 2) MeV. This conclusion should hardly come as a surprise. After all, 
the predominant role played by the model parameters k and A in softening the incompressibility coefficient and g p 
in stiffening the symmetry energy have been known for a long time. What is relevant from the present statistical 
analysis is that we have established quantitatively that the linear Walecka model fails because its prediction for K 
differs from the experimental value by more than 10 standard deviations. We must then conclude that whereas the 
linear Walecka is (fairly) unique, it is not predictive. We now proceed to discuss a particular extension of the Walecka 
model that is highly predictive but not unique: the non-linear FSUGold model. 



B. Example 2: Non-linear FSUGold Model 



Modern relativistic models of nuclear structure have evolved significantly since the early days of the linear Walecka 
model. In the present example we focus on the FSUGold parameter set |18j that is defined by an interacting Lagrangian 
density of the following form: 



9s4>~ 



(<7v^+f r ' tv+la+rs)^) f] $ - ^-*&+^(w»W»y +k(w v W v ) (B M • B") . (29) 

Modifications to the linear Walecka model are motivated by the availability of an ever increasing database of high- 
quality data. For example, the two non-linear scalar terms K and A induce a significant softening of the compression 
modulus of nuclear matter relative to the original Walecka model E] . This is demanded by measurement of 
the giant monopole resonance in medium to heavy nuclei [27]. Further, omega-meson self-interactions, as described 
by the parameter £, also serve to soften the equation of state of symmetric nuclear matter but at much higher 
densities. Indeed, by tuning the value of £ it is possible to produce maximum neutron star masses that differ by 
almost one solar mass while maintaining the saturation properties of nuclear matter intact [9]. Such a softening 
appears consistent with the dynamics of high-density matter as probed by energetic heavy- ion collisions |28j . Finally, 
A v induces isoscalar-isovector mixing and is responsible for modifying the poorly-constrained density dependence of 
the symmetry energy [TTl [^H] • In particular, a softening of the symmetry energy induced by A v appears consistent 
with the distribution of both isoscalar monopole and isovector dipole strength in medium to heavy nuclei [HI [30j [31] . 
In summary, FSUGold is a fairly successful RMF model that has been validated against theoretical, experimental, and 
observational constraints [15]. Note that as additional laboratory and observational data become available (notably 
the recent report of a 2-solar mass neutron star |32j ) refinements to the model my be required |20) . For now, however, 
we will be content with using the FSUGold model to study the small oscillations around the minimum. 

As mentioned earlier, a model should be understood as a combination of an interacting Lagrangian density and a 
quality measure. We define the x 2 " m easure in terms of the following set of observables generated directly from the 
FSUGold parameter set: 



po = 0.1484 fm -3 , 


(30a) 


£ = -16.30 MeV , 


(30b) 


e(2p ) = -5.887 MeV , 


(30c) 


K = 230.0 MeV , 


(30d) 


M* = 0.6100 M , 


(30e) 


J = 26.00 MeV , 


(30f) 


L = 60.52 MeV , 


(30g) 


M max = 1.722 M Q . 


(30h) 



Note that in all cases a 2% uncertainty is attached to all observables — except in the case of the slope of the symmetry 
energy L where the significant larger value of 20% is assumed. This reflects our poor understanding of the density 
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dependence of the symmetry energy. Also note that J represents the value of the symmetry energy at a sub-saturation 
density of psaO.lfm -3 — a density at which the theoretical uncertainties are minimized [33] . Finally, notwithstanding 
the Demorest et al. result G2L the maximum neutron star mass is fixed at M max = 1.722 Mq. Given that a theoretical 



density of psaO.lfm 3 — a density at which the theoretical uncertainties are minimized |33j . Finally, notwithstanding 

tl 

model is used to generate the various observables, a much larger database could be used to define the x 2_m easure, 
if desired. By construction, a very small value for the x 2 -measure is obtained at the FSUGold minimum. We now 
proceed to explore the wealth of information available as one studies deviations around this minimum value. As in 

the previous section, the symmetric matrix of second derivatives M. (now a 7 x 7 matrix) may be diagonalized by 

means of an orthogonal transformation. The diagonal matrix of eigenvalues T> and the matrix of eigenvectors A are 
given by 

V = diag(1.2826x 10 6 , 1.5305 x 10 4 , 4.2472 x 10 2 , 3.2113 x 10 2 , 1.2692 x 10 2 , 6.9619, 3.7690) , (31a) 

/ -7.4967X10" 1 -2.3685X10" 1 3.0853X10" 1 1.2931X10" 1 -5.1254xl0 _1 -8.5089xl0" 2 6.8417xl0" 3 \ 

6.5682X10" 1 -1.5751X10" 1 3.6654X10" 1 1.6504X10" 1 -6.0281X10" 1 -1.3685X lO^ 1 8.5353xl0~ 3 

1.5331 xlO -4 3.1315X10- 3 -7.0050X 10" 1 6.8701xl0~ 2 -3.9206X10' 1 -3.3843xl0~ 2 5.9137X10" 1 

3.8535xl0~ 2 -2.8770X10- 1 -2.4254xl0~ 2 4.7416X10- 1 4.3796xl0" 2 8.2968X10- 1 -5.7643xl0~ 3 

2 _s awiivin-i _i W70vin-i i S77fivin-l i ««vin-i _^ si7fivin-i ofi/iUvm-S 



A- 



3.9417x10"^ -6.8525X10- 1 -1.3772xl0 _i 3.8776X10"- 1 3.5428X10" 1 -4.8376X10" 1 2.6431x10" 
-5.9458X10" 2 6.0558X10" 1 5.4897xl0" 2 7.5689X1Q- 1 6.8021 xlO -2 -2.2175X 10" 1 6.2795xl0" 3 



(31b) 



V 1.2995xl0~ 4 -3.1465xl0~ 3 5.0714X10" 1 -5.7010xl0~ 2 2.9691X10" 1 3.6238xl0~ 2 8.0628X10" 1 
Note that the scaled parameters of the model are associated to the original coupling constants as follows: 

{xi,X2,x 3 ,X4 ) x 5 ,xe,x 7 } -> {g 2 ,g 2 ,g 2 p ,K,\,(,K} ■ (32) 

We observe that the stiffest direction is dominated by two isoscalar parameters and represents — as in the case of the 
linear Walecka model — an out-of-phase oscillation between the scalar attraction and the vector repulsion. Given that 
in RMF models the cancellation between the scalar attraction and the vector repulsion is so delicate, any out-of-phase 
motion yields a significant change in the binding energy per nucleon and a correspondingly dramatic increase in the 
quality measure. The second stiffest direction also involves exclusively isoscalar parameters and is dominated by the 
quartic scalar (A) and vector (£) couplings — and to a lesser extent by the cubic term (k). This linear combination of 
parameters is largely constrained by the incompressibility coefficient Kq and the maximum neutron-star mass M max . 
Although the determination of the maximum neutron-star mass to a 2% accuracy presents a significant observational 
challenge, our statistical analysis suggests that such a determination would strongly constrain the equation of state 
from saturation density up to neutron-star densities. The third stiffest direction (with still a fairly large eigenvalue of 
A3 « 425) is dominated by the two isovector parameters g 2 and A v . For this particular "mode" both parameters oscillate 
out of phase. This behavior can be readily understood by recalling the expression for the symmetry energy |29j : 

E S ym( P ) = + (< - m l + 2A v«) ■ (33) 



P 



In order for the symmetry energy J to remain fixed, then both g 2 and A v must move in phase. If they move out of 
phase, then the symmetry energy can not be kept at this value and the quality measure deteriorates. By the same 
token, the in-phase motion of g 2 p and A v is very poorly constrained — as evinced by the last and softest direction. 
And it is only because the slope of the symmetry energy L was assumed to be somehow constrained (at the 20% 
level) that a positive eigenvalue was even obtained. Note that one of the main goals of the successfully commissioned 
Lead Radius experiment (PREx) at the Jefferson Laboratory is to constrain the density dependence of the symmetry 
energy (i.e., L) by accurately measuring the neutron radius of 208 Pb [331 |3S]. The next to last eigenvalue (Ae ~ 7) 
is also relatively small. This suggest that the out-of-phase motion of the two non-linear scalar couplings (k and A) 
is poorly constrained by the nuclear-matter observables defining the quality measure. Perhaps supplementing the 
quality measure with finite-nuclei observables will help ameliorate this problem. Work along these lines is currently 
in progress. 

We now proceed to estimate theoretical uncertainties as well as to compute correlation coefficients for both the 
model parameters and the physical observables. We start by computing theoretical uncertainties (i.e., variances) for 



the model parameters. These are given by [see Eq. ( 20 )] 



0} = (M- 1 ) = (AV-'A T ) = 5>? A7i (34) 
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and result in the following theoretical uncertainties for the model parameters: 





= 112.19955 ±6.54468 [5.833%] , 


(35a) 


fly 


= 204.54694 ± 15.81183 [7.730%] , 


(35b) 


ft 


= 138.47011 ±42.75427 [30.876%] , 


(35c) 


K 


= 1.42033 ± 0.44827 [31.561%] , 


(35d) 


A 


= 0.02376 ± 0.00445 [18.748%] , 


(35e) 


c 


= 0.06000 ± 0.0057 [9.447%] , 


(35f) 


A v 


= 0.03000 ± 0.01251 [41.711%] . 


(35g) 



We observe that three out of the five isoscalar parameters, namely, g 2 , g 2 , and £, are relatively well constrained (at 
the < 10% level). Whereas g 2 and g 2 are well determined by the saturation properties of symmetric nuclear matter, 
it is the maximum neutron-star mass that constrains C- Yet the remaining two isoscalar parameters (n and A) are 
poorly determined. This is particularly true in the case of k which displays a large (~ 30%) uncertainty. As alluded 
earlier, these large uncertainties develop because the out-of-phase motion of k and A — as controlled by the relatively 
soft sixth eigenvector — is poorly constrained. Given that the in-phase motion of the two isovector parameters (g 2 p and 
A v ) is controlled by the softest of eigenvectors, the theoretical uncertainties in these parameters is also fairly large 
(~30% and «40%, respectively). However, whereas the reason for the latter is associated with the large error bars 
assigned to L, we are unaware at this time on how to better constrain k and A. Perhaps supplementing the quality 
measure with information on various nuclear compressional modes may help resolve this issue. Plans to do so in the 
near future are under consideration. 



CM (O oj > oj o_ 

en en cn y c< jj> 




FIG. 2: (Color online) Color-coded plot of the 21 independent correlation coefficients between the 7 model parameters of the 
FSUGold effective interaction. 

We have computed correlation coefficients between all 21 distinct pairs of model parameters and have displayed 
them in graphical (color-coded) form in Fig. [2] As depicted in the figure, the strongest correlations are between g 2 
and g 2 (0.988), g 2 and A v (0.967), and n and A (-0.962). As alluded in the case of the simpler linear Walecka model, 
the correlations are dominated by the softest directions; in this case the sixth and seventh eigenvectors. Given that 
for these two eigenvectors g 2 and g 2 as well as g 2 and A v move in phase whereas k and A move out of phase, the 
observed correlations ensue. In other words, the three largest eigenvalues strongly suppress the generation of model 
parameters with g 2 and g 2 moving out of phase, k and A in phase, and g 2 and A v out of phase, respectively. Note 
then, that the distribution of isovector parameters g 2 and A v is generated in such a way that the symmetry energy at 
sub-saturation density J remains fixed at 26 MeV (at least within a 2% uncertainty). This quantitative fact validates 
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our heuristic approach — already employed numerous times — to correlate isovector observables (see Refs. [5] [35J 137] 
and references therein). 

We now extend the covariance analysis to the case of physical observables. To do so, we must supply the relevant 
"matrix" of first derivatives [see Eq. (20)]. For completeness we list the first derivatives in tabular form in Table III 
Note that in the case of the model parameters the corresponding matrix of first derivatives is the matrix of eigenvectors 
A. The derivatives encapsulate the sensitivity of the various observables to changes along the different eigenvectors. 
For example, whereas isoscalar observables (such as Eq, pq^ Kq) are insensitive to changes along the mostly isovector 
seventh eigenvector, both L and the neutron-skin thickness of 208 Pb, R n — R p , display a fairly large sensitivity. 
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96 


96 


96 


96 


96 


96 


9e 


1.0551 xl0 +1 


-7.1882X10" 1 


-2.9433 xlO" 2 


4.9029 xlO" 2 


5.1626xl0" 2 


-1.7025 xlO" 2 


-7.7009 xlO" 4 


dpo 


-2.2472 


1.3904 


-9.8868 xlO" 4 


2.2220 xlO" 1 


5.3440 xlO" 2 


1.8998 xlO" 2 


1.5875 xlO- 1 


dK 


-7.4792 


1.3890 


-3.6799 xlO" 2 


-1.8808 xlO" 1 


4.8215xl0" 2 


-2.1691 xlO" 2 


-3.0813 xlO" 4 


DM* 


1.1505 


-5.5581X10" 1 


-1.0614X10" 1 


-8.3170 xlO" 2 


1.6492 xlO" 1 


1.2614xl0" 2 


-2.7701 xlO" 15 


dJ 


-2.7862 xlO" 1 


6.5586xl0" 2 


-3.9136X10" 1 


2.5698 xlO" 2 


-6.8299 xlO" 2 


6.0862 xlO" 4 


6.2323 xlO" 4 


dJ 


-1.6811 


9.4897xl0 _1 


-3.2446 xlO" 1 


1.7826 xlO" 1 


-3.6792 xlO" 2 


6.2151 xlO" 13 


-8.4827xl0" 2 


dL 


-1.8759 


1.1812 


1.5849 xlO" 2 


2.6611 xlO" 1 


-1.0180X10" 1 


-2.2526 xlO" 2 


-3.8593 xlO" 1 


d(Rn-R P ) 


7.0224 


9.6188xl0" 2 


-2.3362 xlO -1 


1.9804 xlO" 1 


-2.5651 xlO" 2 


-1.9368 xlO" 2 


-3.4167X10" 1 


dRi.o 


9.9947 xlO" 1 


-3.1989X10" 1 


-1.6534 xlO" 2 


-6.5939 xlO" 2 


-4.5197xl0" 2 


-4.5964 xlO" 13 


-3.9800 xlO" 2 


dRi.4 


5.0300 xlO" 1 


-3.0884xl0 _i 


9.3778 xlO" 3 


-1.1119X10" 1 


-6.3792 xlO" 2 


3.1016 x 10"* 3 


-2.6571 xlO" 2 


9M max 


-2.7675 xlO" 1 


-1.4882X10" 1 


3.1173xl0" 2 


-1.6394X10" 1 


-6.8790 xlO" 2 


3.4817xl0" 2 


-2.4367 xlO" 15 



TABLE III: First derivatives of the scaled observables (i.e., observable scaled to its value at the x 2 -minimum) as a function of 
6 at the x 2 -minimum; see Eq. (20 1. 



Given the enormous interest in constraining the density dependence of the symmetry energy, we estimate theoretical 
uncertainties on three — mostly isovector — observables. These are the symmetry energy at saturation density J, the 
neutron-skin thickness of 208 Pb, and the radius (Ria) of a M = 1.4M Q neutron star. Recall that it was J (not J) 
that was included in the definition of the quality measure. We obtain, 

J = (32.593 ± 1.574) MeV [4.830%] , (36a) 
R n -Rp = (0.207 ± 0.037) fm [17.698%] , (36b) 
Ria = (11.890 ± 0.194) km [1.631%] . (36c) 

We now comment on each of these cases individually. Before we do so, however, note that correlation coefficients for 
11 observables (i.e., 55 independent pairs) are depicted in a color-coded format in Fig. [3] First, the central value of J 
along with its theoretical uncertainty may be easily understood by invoking a first-order expansion for the symmetry 
energy J at sub-saturation density (/5o« 0.103 fm" 3 ) in terms of J and L [38]. That is, 

J = J + xL + ...P3 (32.208 ± 1.346) MeV , x = 1 [ 1 - ^) « 0.103 , (37) 

3 V PoJ 

where the errors were added in quadrature. So although J is strongly correlated to L (with a correlation coefficient 
of 0.922) the error in the former is significantly smaller than the latter because of the small value of x. Second, for 
the neutron-skin thickness of 208 Pb we find a theoretical error comparable to the one assumed for L and a correlation 
coefficient between the two observables of almost one (0.995). Such a strong correlation is consistent with two recent 
studies that employ a large number of accurately-calibrated relativistic and non-relativistic interactions to uncover 
the correlation [39l [40] . Also consistent with these studies, specifically with Ref. [40], is the fact that the proposed 1% 
measurement of the neutron radius of 208 Pb by the PREx collaboration [331 [35] may not be able to place a significant 
constrain on L. For example, our covariance analysis suggests that the 20% uncertainty assumed for L translates 
into a theoretical error in the neutron skin of 0.037 fm — or about a 0.7% uncertainty in the neutron radius of 208 Pb. 
Conversely, if L is to be determined to within 10% (i.e., L«60 ± 6 MeV) then the neutron skin must be constrained 
to about 0.018 fm so the neutron radius must be measured with close to a 0.3% accuracy — a fairly daunting task. 
Finally, we obtain a very small theoretical uncertainty for the radius of a 1.4 solar-mass neutron star and a correlation 
coefficient between L and R\a (or R n — R p and Ria) of 0.811. Although the radius of the neutron star is sensitive 
to the density dependence of the symmetry energy }41j . Ria can not be uniquely constrained by a measurement of 
R n — R p because whereas the latter depends on the symmetry energy at (or below) saturation density, the former 
is also sensitive to the equation of state at higher densities [2H]- Note that a far better correlation coefficient of 
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FIG. 3: (Color online) Color-coded plot of the 55 independent correlation coefficients between 11 physical observables as 
computed with the FSUGold effective interaction. 

0.942 is obtained between L and the radius of a 1.0 solar-mass neutron star. Regardless (with all things being equal) 
knowledge of the slope of the symmetry energy to a 20% accuracy significantly constrains the stellar radius. 

IV. CONCLUSIONS 

The demand for theoretical predictions that include meaningful and reliable uncertainties is increasing. Such a 
sentiment has been articulated in a recent publication by the editors of the Physical Review A [2]. The need to 
quantify model uncertainties in an area such as theoretical nuclear physics is particularly urgent as models that are 
fitted to experimental data are then used to extrapolate to the extremes of temperature, density, isospin asymmetry, 
and angular momentum. Inspired by some of the central ideas developed in Ref. [5J, a systematic statistical approach 
was applied to a class of relativistic mean-field models. The aim of this statistical analysis was twofold. First, to 
attach meaningful and reliable theoretical uncertainties to both the model parameters as well as to the predicted 
observables. Second, to quantify the degree of correlation between physical observables. 

Modern relativistic mean-field models have evolved considerably since the early days of the linear Walecka model. 
Based on certain shortcoming of the Walecka model — most notably the inability to reproduce the incompressibility 
coefficient of symmetric nuclear matter — the Lagrangian density was augmented by non-linear cubic and quartic 
scalar-meson terms. However, based on modern effective-field-theory tenets, such as naturalness and power counting, 
a consistent Lagrangian density should include all terms up to fourth order in the meson fields. But in doing so, how 
should one constrain the large number of model parameters? In principle, one should follow the standard protocol 
of determining all model parameters through a ^-minimization procedure. In practice, however, many successful 
theoretical approaches arbitrarily set some of the model parameters to zero. The argument behind this fairly ad- 
hoc procedure is that the full set of parameters is poorly determined by existing data, so ignoring a subset model 
parameters does not compromise the quality of the fit. 

A covariance analysis such as the one implemented here should be able to clarify in a quantitative fashion the precise 
meaning of a "poorly determined set of parameters" . To do so, one should focus — not on the minimum of the x 2 ~ 
measure but rather — on its behavior around the minimum. As in any small-oscillations problem, the deviations around 
the minimum are controlled by a symmetric matrix of second derivatives that may be used to extract theoretical error 
bars and to compute correlation coefficients among physical observables. However, to access the wealth of information 
available in the covariance analysis we took it a step further and diagonalized the matrix of second derivatives. Upon 
diagonalization, the deviations of the x 2 -measure from the minimum are parametrized in terms of a collection of 
"uncoupled harmonic oscillators" . By doing so, one could readily identify stiff and soft modes in parameter space, 
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namely, eigenvectors characterized by either large or small eigenvalues, respectively. 

We now summarize some of the most important lessons learned. First, a stiff direction represents a particular 
linear combination of model parameters that is well constrained by the set of physical observables included in the 
X 2 -measure. By the same token, a soft direction suggests that additional physical observables are required to further 
constrain the model. Second, given that model parameters around the minimum are distributed according to the 
X 2 -measure, the soft directions dominate the correlation analysis. Finally, testing whether a model is well constrained 
by individually varying its parameters — rather than by varying them coherently as suggested by the structure of 
the eigenvectors — may be misleading. To illustrate these findings we used two relatively simple, yet illuminating, 
examples: (a) the linear Walecka model and (b) the FSUGold parametrization. Note that ultimately we aim to 
implement the covariance analysis with a x 2 -measure defined by a consistent Lagrangian density. 

A particularly clear example of a stiff direction was represented by the out-of-phase motion of the scalar g s and 
vector g v coupling constants in the linear Walecka model. Indeed, increasing the scalar attraction while at the same 
time reducing the vector repulsion leads to a significant increase in the binding energy per nucleon and, thus, in a 
significant deterioration of the x 2 -measure. The in-phase motion of g s and g v , however, is not as well constrained (the 
ratio of the two eigenvalues is about 1000). Therefore, configurations in parameter space generated by the x 2 -measure 
were dominated by pairs of coupling constants that were in phase, thereby resulting in a correlation coefficient between 
g s and g v that was, as expected, large and positive. Note, however, that if g s and g v were varied individually, one would 
erroneously conclude that the model is much better constrained than it really is — since changes in the x 2_m easure 
would be dominated by the largest eigenvalue. 

In our second example we considered the accurately-calibrated FSUGold interaction with an isovector interaction 
determined by two parameters (g p and A v ). We found the out-of-phase motion of g p and A v to be strongly constrained 
by the value of the symmetry energy at a density of about 0.1 fm. However, our poor knowledge of the density 
dependence of the symmetry energy left the in-phase motion of g v and A v largely unconstrained. Effectively then, 
correlations in the isovector sector were induced by the in-phase motion of g p and A v — subject to the constraint that 
the symmetry energy at p«0.1 fm remains intact. This procedure validates the heuristic approach that we have used 
for some time to estimate correlations among isovector observables. Yet a benefit of the present analysis is that one 
can precisely quantify the theoretical errors as well as the correlation among observables. For example, we concluded 
that if the slope of the symmetry energy is to be determined with a 10% uncertainty, then the neutron-skin thickness 
of 208 Pb should be measured with a 0.3% accuracy. This more stringent limit seems to agree with the conclusions of 
Ref. 0D|. 

In the future we aim to apply the covariance analysis discussed here to the construction of a relativistic density 
functional that will include all terms up to fourth order in the meson fields. Moreover, we plan to calibrate the x 2 - 
measure using various properties of finite nuclei and neutron stars. In addition, we reiterate a point made in Ref. [5] 
that the methodology used in this work should be applicable to any problem where model parameters are determined 
from optimizing a quality measure. 
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